################################################################################
# CONGRESSIONAL INTERVENTION IN AGENCY ADJUDICATION
# AUTHORS: LINDSEY GAILMARD, DANIEL E. HO, MARK KRASS
################################################################################

################################################################################
# THIS FILE CREATES TABLES AND FIGURES PRESENTED IN THE MAIN PAPER:
# FIGURE 3: CONGRESSIONAL INTERVENTION BY YEAR
# TABLE 1: BALANCE TABLE
# TABLE 2: AVERAGE DISTRICT-LEVEL DEMOGRAPHIC CHARACTERISTICS
# FIGURE 4: CONGRESSIONAL INTERVENTION RATE BY DISTRICT
# FIGURE 5: NUMBER OF INTERVENTIONS BY IAVA GRADE
# FIGURE 6: SHARE OF APPEALS IN STUDY WITH ADVANCEMENT BY DECISION YEAR
# FIGURE 7: DISTRIBUTION OF APPEAL DURATION BY CONGRESSIONAL INQUIRY STATUS
# FIGURE 8: ASMD FOR MATCHED SAMPLE
# FIGURE 9: FRACTION OF APPEALS ADVANCED ON DOCKET WITHOUT DOCUMENTED CAUSE
# FIGURE 10: KAPLAN-MEIER CURVE FOR ADJUSTED TIME TO FINAL DECISION
# FIGURE 12: MEDIATION ANALYSIS
################################################################################

# LOAD PACKAGES
# install.packages("tidyverse")
# install.packages("extrafont")
# install.packages("readxl")
# install.packages("svMisc")
# install.packages("GenMatch")
# install.packages("Matching")
# install.packages("magick")
# install.packages("cem")
# install.packages("rstatix")
install.packages("DiagrammeR")
install.packages("fixest")
install.packages("modelsummary")
install.packages("readr")

library(DiagrammeR)
library(mediation)
library(fixest)
library(modelsummary)
library(rstatix)
library(readxl)
library(readr)
library(extrafont)
require(tigris)
require(ggplot2)
require(sf)
require(ggmap)
require(foreign)
require(RCurl)
require(xlsx)
require(zipcode)
require(choroplethr)
require(dplyr)
require(tidycensus)
library(lubridate)
require(RColorBrewer)
require(readr)
require(maps)
require(maptools)
require(data.table)
require(kableExtra)
require(stringr)
require(tidyverse)
require(gridExtra)
require(survival)
require(survminer)
require(GenMatch)
require(Matching)
require(MatchIt)
require(svMisc)
require(cem)

library(tidyverse)

font_family <- 'sans'

tryCatch({
  font_import(pattern = 'Times New Roman', prompt = FALSE)
  font_family <- 'Times New Roman'
})

# SET LOCAL WORKING DIRECTORY
setwd("/Users/lindseygailmard/Dropbox/RegLab/BVA Congressional Intervention YLJ")

getwd()

# LOAD ALL DATASETS FOR FIGURES AND TABLES
fig3 <-
  read.csv("./Replication Files/Data/fig3.csv", header = TRUE)
table1 <-
  read.csv("./Replication Files/Data/table1.csv", header = TRUE)
table2_fig4 <-
  read.csv("./Replication Files/Data/table2 - fig4.csv", header = TRUE)
iava_110 <-
  read.csv("./Replication Files/Data/fig5_a.csv", header = TRUE)
iava_111 <-
  read.csv("./Replication Files/Data/fig5_b.csv", header = TRUE)
matched_dataset <-
  read.csv("./Replication Files/Data/matched dataset.csv", header = TRUE)

################################################################################
############## FIGURE 3: CONGRESSIONAL INTERVENTION BY YEAR  ###################
################################################################################

ggplot(data = fig3, aes(x = year, y = rate)) + geom_line() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    text = element_text(size = 14),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    panel.border = element_rect(
      colour = "black",
      fill = NA,
      size = 0.75
    ),
    axis.line = element_line(colour = "black")
  ) +
  ylim(0, 0.15) + xlab("Year") + ylab("Congressional Intervention/New Appeals") +
  ggtitle("Congressional Intervention by Year (Rate)")  +
  scale_y_continuous(breaks = c(0.05, 0.10, 0.15),
                     sec.axis = sec_axis( ~ . * 1, breaks = c(0.05, 0.10, 0.15))) +
  annotate(
    "rect",
    xmin = 2003,
    xmax = 2017,
    ymin = -Inf,
    ymax = Inf,
    alpha = 0.25
  ) +
  annotate("text",
           x = 2010,
           y = 0.15,
           label = "Study Period")

ggsave(
  "./Replication Files/Figures/Figure 3.png",
  device = "png",
  width = 5,
  height = 4
)

################################################################################
##########################  TABLE 1: BALANCE TABLE #############################
################################################################################

# T-TESTS FOR DIFFERENCES BETWEEN APPEALS WITHOUT/WITH CI
# APPLIES BENJAMINI HOCHBERG CORRECTION TO P-VALUES FOR MULTIPLE TESTS
runTTests <- function(data, outcome.var.names, indp.var) {
  #############################
  # data = DATAFRAME CONTAINING COLUMNS FOR split.pop.var, outcome.var, and indp.var
  # outcome.var.names = CHARACTER VECTOR OF VARIABLE NAMES USED FOR OUTCOMES
  # indp.var = STRING NAME OF VARIABLE USED FOR INDEPENDENT VARIABLE
  #############################
  # PREP LOOP
  results.df <- c()
  # SET COLUMN NAMES
  col.names <- c("outcome",
                 "control_rate",
                 "ci_rate",
                 "t_value",
                 "p_value",
                 "df",
                 "se.diff")
  for (outcome.var in outcome.var.names) {
    # T-TEST FOR CONTINUOUS VARIABLES
    possibleError <- tryCatch(
      ttest <- t.test(data[[outcome.var]] ~ data[[indp.var]]),
      error = function(e)
        e
    )
    if (inherits(possibleError, "error")) {
      tcount <- table(data[[indp.var]])
      test.row <- c(outcome.var, NA, NA, NA, NA, NA, NA)
    } else{
      tcount <- table(data[[indp.var]])
      test.row <- c(
        outcome.var,
        ttest$estimate[1],
        ttest$estimate[2],
        ttest$statistic,
        ttest$p.value,
        ttest$parameter[['df']],
        (ttest$estimate[1] - ttest$estimate[2]) / ttest$statistic
      )
    }
    # SET NAMES
    names(test.row) <- col.names
    results.df <- rbind(results.df, test.row)
  }
  results.df <-
    suppressWarnings(data.frame(results.df, stringsAsFactors = F))
  names(results.df) <- c("outcome",
                         "control_rate",
                         "ci_rate",
                         "t_value",
                         "p_value",
                         "df",
                         "se.diff")
  # CONVERT TO NUMERIC
  results.df[, c(2:ncol(results.df))] <-
    sapply(results.df[, c(2:ncol(results.df))], as.numeric)
  results.df[, 2:5] <- round(results.df[, 2:5], 2)
  results.df[, 7] <- round(results.df[, 7], 2)
  
  return(results.df)
}

outcomes = c(
  "male",
  "appellant_age_nod",
  "n_issues",
  "finhard",
  "termill",
  "advage",
  "no_aod_criteria",
  "ctyp_ps3",
  "ctyp_ps4",
  "ctyp_ps5",
  "ctyp_ps6",
  "ctyp_ps7",
  "ctyp_ps8",
  "ctyp_ps9",
  "rep_unrepresented",
  "rep_attorney",
  "rep_servorg",
  "mil",
  "female",
  "house_incumbent",
  "rep",
  "legal",
  "VA_committee",
  "ways_and_means_committee",
  "appropriations_committee",
  "budget_committee",
  "medianincome",
  "medianage",
  "prcntunemp",
  "prcntnotemploy",
  "prcnths",
  "prcntba",
  "prcntblack",
  "prcntwhite"
)

# RUN BALANCE TEST
balance.tests <- runTTests(data = table1,
                           outcome.var.names = outcomes,
                           indp.var = "CI")
# ADJUST FOR MULTIPLE TESTS
balance.tests$p_value_adjusted <- NA
balance.tests$p_value_adjusted <-
  p.adjust(as.numeric(balance.tests$p_value), method =
             "fdr")
balance.tests$p_value_adjusted <-
  round(balance.tests$p_value_adjusted, 2)

rownames(balance.tests) <- c(
  "Appellant is male",
  "Appellant Age at NOD (Years)",
  "Issues Per Appeal",
  "Financial Hardship",
  "Terminal Illness",
  "Advanced Age",
  "No AOD Criteria Documented",
  "WWII (9/16/40-7/25/47)",
  "Peacetime (7/26/47-6/26/50)",
  "Korean Conflict (6/27/50-1/31/55)",
  "Post-Korea (2/1/55-8/4/64)",
  "Vietnam Era (8/5/64- 5/7/75)",
  "Post-Vietnam (5/8/75-8/1/90)",
  "Persian Gulf (8/2/90-Present)",
  "Unrepresented",
  "Attorney",
  "Service Org or Agency",
  "Veteran",
  "Female",
  "Incumbent",
  "Republican",
  "Legal Background",
  "Veterans' Affairs",
  "Ways and Means",
  "Appropriations",
  "Budget",
  "Median Income",
  "Median Age",
  "Unemployment Rate (%)",
  "Population Out of the Labor Force (%)",
  "High School Education (%)",
  "College Education (%)",
  "Black Population",
  "White Population"
)
balance.tests <- balance.tests[, c(2:3, 5, 8)]
colnames(balance.tests) <- c("No Intervention",
                             "Congressional Intervention",
                             "p-value",
                             "adjusted p-value")

# CREATE TABLE 1
kbl(balance.tests, format = "latex", booktabs = T) %>%
  kable_styling(latex_options = "scale_down") %>%
  pack_rows("Appellant Demographic Characteristics", 1, 2) %>%
  pack_rows("Advancement on the Docket Criteria", 3, 6) %>%
  pack_rows("Period of Service", 7, 13) %>%
  pack_rows("Veteran Representation", 14, 17) %>%
  pack_rows("House Representative Characteristics", 18, 22) %>%
  pack_rows("House Committee Membership", 23, 26) %>%
  pack_rows("Congressional District Characteristics", 27, 34)

# TOTAL NUMBER OF UNIQUE APPELLANTS
length(unique(table1[table1$CI == 0, 'appellant_id']))
length(unique(table1[table1$CI == 1, 'appellant_id']))

# TOTAL NUMBER OF UNIQUE APPEALS
length(unique(table1[table1$CI == 1, 'appeal_id']))
length(unique(table1[table1$CI == 0, 'appeal_id']))

################################################################################
########### FIGURE 4: CONGRESSIONAL INTERVENTION RATE BY DISTRICT ##############
################################################################################

us_states <- states(cb = TRUE, resolution = "20m") %>%
  filter(NAME != "Puerto Rico") %>%
  shift_geometry()

us_states <- us_states %>%
  rename(state_fips = STATEFP)

# CONGRESSIONAL DISTRICT SHAPE FILES
# https://cdmaps.polisci.ucla.edu/

get_congress_map <- function(cong = 113) {
  tmp_file <- tempfile()
  tmp_dir  <- tempdir()
  zp <- sprintf("http://cdmaps.polisci.ucla.edu/shp/districts%03i.zip",
                cong)
  download.file(zp, tmp_file)
  unzip(zipfile = tmp_file, exdir = tmp_dir)
  fpath <- paste(tmp_dir,
                 sprintf("districtShapes/districts%03i.shp", cong),
                 sep = "/")
  st_read(fpath)
}

cd114 <- get_congress_map(114)

cd114 <- merge(cd114, fips_codes_state, by = "STATENAME")
cd114$DISTRICT <- as.numeric(cd114$DISTRICT)
cd114$cd_match <- cd114$state_fips * 100 + cd114$DISTRICT

cd114 <- cd114 %>%
  shift_geometry()

intervention_114 <- subset(table2_fig4, congress == 114)

cd114_intervention <- cd114 %>% left_join(intervention_114, by = "cd_match")
cd114_intervention$intervention <- ifelse(is.na(cd114_intervention$district_CI),
                                          0,
                                          cd114_intervention$district_CI)
cd114_intervention$intervention_pct <- 100 * (cd114_intervention$intervention /
                                                cd114_intervention$Veteran.Population)

ggplot() +
  geom_sf(
    data = cd114_intervention,
    aes(fill = intervention_pct),
    inherit.aes = FALSE,
    alpha = 0.9,
    lwd = 0
  ) +
  geom_sf(data = us_states, alpha = 0) +
  scale_fill_distiller(palette = "Blues", direction = 1) +
  theme_void() +
  labs(fill = 'Intervention Rate (%)') +
  theme_void() +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    text = element_text(family = "Times New Roman", size = 12),
    plot.caption.position = "plot",
    plot.caption = element_text(
      hjust = 0,
      family = "Times New Roman",
      size = 12
    )
  )

ggsave(
  "./Replication Files/Figures/Figure 4.png",
  device = "png",
  width = 6.5,
  height = 3
)

################################################################################
#########  TABLE 2: AVERAGE DISTRICT-LEVEL DEMOGRAPHIC CHARACTERISTICS #########
################################################################################

table2 <- table2_fig4 %>%
  filter(congress == 110 | congress == 112) %>%
  mutate(intervention_pct = 100 * district_CI / Veteran.Population) %>%
  rename(intervention = district_CI)

table2$intervention_pct <- ifelse(is.na(table2$intervention_pct), 0, table2$intervention_pct)

table2$intervention <- ifelse(is.na(table2$intervention), 0, table2$intervention)

summary(table2$intervention[table2$congress == 110])
summary(table2$intervention[table2$congress == 112])

pct75_110 <- table2 %>% filter(congress == 110) %>% pull(intervention) %>% quantile(probs = .75)
pct75_112 <- table2 %>% filter(congress == 112) %>% pull(intervention) %>% quantile(probs = .75)

table2$top <- ifelse(
  table2$intervention >= pct75_110 & table2$congress == 110,
  1,
  ifelse(table2$intervention >= pct75_112 & table2$congress == 112, 1, 0)
)

table2_summary <- table2 %>%
  group_by(congress, top) %>%
  summarise(
    income = mean(medianincome),
    age = mean(medianage),
    unemp = mean(prcntunemp),
    notemploy = mean(prcntnotemploy),
    ba = mean(prcntba),
    hs = mean(prcnths),
    black = mean(prcntblack),
    white = mean(prcntwhite)
  )

################################################################################
############## FIGURE 5: NUMBER OF INTERVENTIONS BY IAVA GRADE #################
################################################################################

iava_111 <- iava_111 %>%
  filter(grepl("^[A-Za-z]", STATENAME.x))
iava_111 <- iava_111 %>% dplyr::select(
  STATENAME.x,
  DISTRICT.x,
  party.y,
  name,
  sort_name,
  chamber,
  IAVA_2010_GRADE,
  IAVA_2010_SCORE,
  district_CI
) %>%
  mutate(assessment_yr = "111th Congress\n(2010 Report)", cong = 111) %>%  # Adds a new column 'assessment_yr' with the value 2010
  rename(
    IAVA_SCORE = IAVA_2010_SCORE,
    IAVA_GRADE = IAVA_2010_GRADE,
    STATENAME = STATENAME.x,
    DISTRICT = DISTRICT.x,
    party = party.y,
    intervention = district_CI
  )

iava_110 <- iava_110 %>%
  filter(grepl("^[A-Za-z]", STATENAME.x))
iava_110 <- iava_110 %>% dplyr::select(
  STATENAME.x,
  DISTRICT.x,
  party.y,
  name,
  sort_name,
  chamber,
  Iava_2008_Grade,
  Iava_2008_Score,
  district_CI
) %>%
  mutate(assessment_yr = "110th Congress\n(2008 Report)", cong = 110) %>%  # Adds a new column 'assessment_yr' with the value 2010
  rename(
    IAVA_SCORE = Iava_2008_Score,
    IAVA_GRADE = Iava_2008_Grade,
    STATENAME = STATENAME.x,
    DISTRICT = DISTRICT.x,
    party = party.y,
    intervention = district_CI
  )

colnames(iava_110) == colnames(iava_111)

iava_110$IAVA_SCORE <- as.numeric(iava_110$IAVA_SCORE)

iava_data <- bind_rows(iava_110, iava_111)

pct75 <- iava_data %>% pull(intervention) %>% quantile(probs = .75, na.rm = T)
pct50 <- iava_data %>% pull(intervention) %>% quantile(probs = .50, na.rm = T)

iava_data %>% filter(iava_data$intervention > pct75) %>%
  ggplot(aes(x = IAVA_SCORE, y = intervention, color = party)) +
  geom_point() +              # Creates the scatterplot
  geom_smooth(method = "lm") +  # Adds a smoothed trendline using LOESS method
  facet_wrap( ~ assessment_yr) +
  labs(x = "111th Congress\n(2008 Report)", y = "Intervention", title = "Scatterplot with Smoothed Trendline") +
  theme_minimal()

data.filtered <- iava_data %>% filter(intervention > pct75)

cor.test(iava_data$IAVA_SCORE, iava_data$intervention, method = "spearman")
cor.test(data.filtered$IAVA_SCORE, data.filtered$intervention)

iava_data$grades <-  factor(iava_data$IAVA_GRADE, levels = c("F", "D", "C", "B", "A", "A+"))
iava_data$two_grades <-  ifelse(iava_data$IAVA_GRADE %in% c("A", "A+"), "A/A+", "All Others")

plt <- iava_data %>% filter(!is.na(two_grades)) %>%
  ggplot(aes(x = two_grades, y = intervention)) +
  geom_boxplot() +              # Creates the boxplot
  facet_wrap( ~ assessment_yr) +
  labs(x = "IAVA Grade", y = "Annual Interventions", title = "") +
  theme(
    text = element_text(family = "serif"),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
  )

ggsave(
  "./Replication Files/Figures/Figure 5.png",
  device = "png",
  width = 4.5,
  height = 3
)

result <- iava_data %>%
  # Create a new column to define the quadrants
  mutate(
    quadrant = case_when(
      cong == 111 & two_grades == "All Others" ~ "111 All Others",
      cong == 111 & two_grades != "All Others" ~ "111 A/A+",
      cong == 110 & two_grades == "All Others" ~ "110 All Others",
      cong == 110 & two_grades != "All Others" ~ "110 A/A+",
      TRUE ~ "Other"  # Handle any other cases if necessary
    )
  ) %>%
  # Group by the new quadrant
  group_by(quadrant) %>%
  # Calculate the median of the `interventions` column for each quadrant
  summarize(median_interventions = median(intervention, na.rm = TRUE))

# Print the result
print(result)

result.party <- dat %>%
  mutate(isA = ifelse(two_grades == "All Others", 0, 1)) %>%
  # Create a new column to define the quadrants
  mutate(quadrant = case_when(
    party == "republican" ~ "rep",
    party != "republican" ~ "dem",
    TRUE ~ "Other"  # Handle any other cases if necessary
  )) %>%
  # Group by the new quadrant
  group_by(quadrant) %>%
  # Calculate the median of the `interventions` column for each quadrant
  summarize(
    median_interventions = median(intervention, na.rm = TRUE),
    iava = mean(isA, na.rm = TRUE)
  )
print(result.party)


################################################################################
##### FIGURE 6: SHARE OF APPEALS IN STUDY WITH ADVANCEMENT BY DECISION YEAR ####
################################################################################

fig6 <- table1 %>%
  filter(decyear > 2002 & decyear < 2018) %>%
  group_by(decyear) %>%
  summarise(mean_aod = mean(aod))

ggplot(data = fig6, aes(x = decyear, y = mean_aod)) +
  geom_bar(stat = "identity",
           width = 0.5,
           position = position_dodge(0.7), 
           fill = "dark gray") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 12),
    text = element_text(size = 12),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = "black")
  ) +
  ylim(0, 0.1) + xlab("Year") + ylab("Share of Appeals Decided with AOD") +
  ggtitle("")

ggsave(
  "./Replication Files/Figures/Figure 6.png",
  device = "png",
  width = 9,
  height = 4
)

################################################################################
### FIGURE 7: DISTRIBUTION OF APPEAL DURATION BY CONGRESSIONAL INQUIRY STATUS ##
################################################################################

# DURATION OF APPEALS
# REMOVE APPEALS WITH NEGATIVE OR 0 FOR DURATION
CI_only <- subset(table1, CI == 1)
nCI_only <- subset(table1, CI == 0)

CI_only <- CI_only[!(is.na(CI_only$bfddec)), ]
CI_only <- CI_only[!(is.na(CI_only$bfdnod)), ]

CI_only <- CI_only %>%
  mutate(duration = interval(noddate, decdate) %>%
           time_length(unit = "days") %>%
           trunc())

CI_only$duration <- as.numeric(CI_only$duration)
CI_only$duration_years <- CI_only$duration / 365
neg <-  subset(CI_only, duration_years <= 0)
less1year <- subset(CI_only, duration_years < 1)

nCI_only <- nCI_only[!(is.na(nCI_only$bfddec)), ]
nCI_only <- nCI_only[!(is.na(nCI_only$bfdnod)), ]

nCI_only <- nCI_only %>%
  mutate(duration = interval(noddate, decdate) %>%
           time_length(unit = "days") %>%
           trunc())

nCI_only$duration <- as.numeric(nCI_only$duration)
nCI_only$duration_years <- nCI_only$duration / 365
neg <- subset(nCI_only, duration_years <= 0)
less1yearnCI <- subset(nCI_only, duration_years < 1)

CI_only <- subset(CI_only, duration_years > 0)
CI_only$duration_years_round <- round(CI_only$duration_years, 0)
nCI_only <- subset(nCI_only, duration_years > 0)
CI_only$duration_years_round <- ifelse(
  CI_only$duration_years < 1,
  1,
  ifelse(
    CI_only$duration_years >= 1 & CI_only$duration_years < 2,
    2,
    ifelse(
      CI_only$duration_years >= 2 & CI_only$duration_years < 3,
      3,
      ifelse(
        CI_only$duration_years >= 3 & CI_only$duration_years < 4,
        4,
        ifelse(
          CI_only$duration_years >= 4 & CI_only$duration_years < 5,
          5,
          ifelse(
            CI_only$duration_years >= 5 & CI_only$duration_years < 6,
            6,
            ifelse(
              CI_only$duration_years >= 6 & CI_only$duration_years < 7,
              7,
              ifelse(
                CI_only$duration_years >= 7 & CI_only$duration_years < 8,
                8,
                ifelse(CI_only$duration_years >= 8 &
                         CI_only$duration_years < 9, 9, 10)
              )
            )
          )
        )
      )
    )
  )
)
NCI_only = nCI_only

NCI_only$duration_years_round <-
  ifelse(
    NCI_only$duration_years < 1,
    1,
    ifelse(
      NCI_only$duration_years >= 1 & NCI_only$duration_years < 2,
      2,
      ifelse(
        NCI_only$duration_years >= 2 & NCI_only$duration_years < 3,
        3,
        ifelse(
          NCI_only$duration_years >= 3 & NCI_only$duration_years < 4,
          4,
          ifelse(
            NCI_only$duration_years >= 4 & NCI_only$duration_years < 5,
            5,
            ifelse(
              NCI_only$duration_years >= 5 & NCI_only$duration_years < 6,
              6,
              ifelse(
                NCI_only$duration_years >= 6 & NCI_only$duration_years < 7,
                7,
                ifelse(
                  NCI_only$duration_years >= 7 & NCI_only$duration_years < 8,
                  8,
                  ifelse(NCI_only$duration_years >= 8 &
                           NCI_only$duration_years < 9, 9, 10)
                )
              )
            )
          )
        )
      )
    )
  )

ci <- round(table(CI_only$duration_years_round) / dim(CI_only)[1], 3)
ci <- t(ci)

nci <-
  round(table(NCI_only$duration_years_round) / dim(NCI_only)[1], 3)
nci <- t(nci)

t <- rbind(ci, nci)
rownames(t) <-
  c("Congressional Intervention", "No Congressional Intervention")
colnames(t) <-
  c("<1", "1-2", "2-3", "3-4", "4-5", "5-6", "6-7", "7-8", "8-9", "10+")

png(
  file = "./Replication Files/Figures/Figure 7.png",
  width = 9,
  height = 4
)

par(
  mar = c(4, 4, 2, 1),
  mgp = c(.7, 0.1, 0),
  tcl = -0.3,
  cex.lab = .8,
  cex.main = 1,
  cex.axis = .6
)
barplot(
  t,
  main = "Distribution of Appeal Duration",
  xlab = "Duration of Appeal (Years)",
  col = c("black", "gray"),
  beside = TRUE,
  ylim = c(0, .3),
  ylab = "Proportion of Appeals"
)
legend(
  "topright",
  legend = rownames(t),
  fill = ,
  pch = 15,
  col = c("black", "gray"),
  ncol = 1,
  cex = 0.8,
  bty = "n"
)

dev.off()

################################################################################
###################### FIGURE 8: ASMD FOR MATCHED SAMPLE #######################
################################################################################

balance_all <- dplyr::select(
  table1,
  c(
    "CI",
    "male",
    "appellant_age_nod",
    "ctyp_ps0",
    "ctyp_ps1",
    "ctyp_ps2",
    "ctyp_ps3",
    "ctyp_ps4",
    "ctyp_ps5",
    "ctyp_ps6",
    "ctyp_ps7",
    "ctyp_ps8",
    "ctyp_ps9",
    "rep_unrepresented",
    "rep_attorney",
    "rep_servorg"
  )
)

# Transform the data into long format
# Put all variables in the same column except `CI`, the grouping variable
balance_data.long <- balance_all %>%
  pivot_longer(-CI, names_to = "variables", values_to = "value")

smd.test <- balance_data.long %>%
  group_by(variables) %>%
  cohens_d(value ~ CI, var.equal = TRUE)

balance_matched <- dplyr::select(
  matched_dataset,
  c(
    "CI",
    "male",
    "appellant_age_nod",
    "ctyp_ps0",
    "ctyp_ps3",
    "ctyp_ps4",
    "ctyp_ps5",
    "ctyp_ps6",
    "ctyp_ps7",
    "ctyp_ps8",
    "ctyp_ps9",
    "rep_unrepresented",
    "rep_attorney",
    "rep_servorg"
  )
)

# Transform the data into long format
# Put all variables in the same column except `CI`, the grouping variable
matched.long <- balance_matched %>%
  pivot_longer(-CI, names_to = "variables", values_to = "value")

smd.test.matched <- matched.long %>%
  group_by(variables) %>%
  cohens_d(value ~ CI, var.equal = TRUE)

# CREATE TABLE OF SMD FOR BALANCE PLOT
smd.balance <- dplyr::select(smd.test.matched, c("variables", "effsize"))
smd.balance$post.match.smd <- abs(smd.balance$effsize)
smd.balance <- dplyr::select(smd.balance, c("variables", "post.match.smd"))

smd.pre <- dplyr::select(smd.test, c("variables", "effsize"))
smd.pre$abs_effsize <- abs(smd.pre$effsize)
smd.pre$sign <- ifelse(smd.pre$effsize > 0, 1, 0)
smd.balance <- left_join(smd.balance, smd.pre, by = "variables")

smd.balance <- smd.balance %>%
  mutate(
    variables = recode(
      variables,
      male = "Male",
      appellant_age_nod = "Age at NOD",
      rep_unrepresented = "Unrepresented",
      rep_attorney = "Attorney",
      rep_servorg = "Service Organization/State Agency",
      ctyp_ps0 = "No Creditable Service",
      ctyp_ps3 = "WWII Service",
      ctyp_ps4 = "Peacetime Service",
      ctyp_ps5 = "Korean Conflict Service",
      ctyp_ps6 = "Post-Korea Service",
      ctyp_ps7 = "Vietnam Era Service",
      ctyp_ps8 = "Post-Vietnam Service",
      ctyp_ps9 = "Persian Gulf Service"
    )
  )

data2 <- smd.balance %>%
  dplyr::select(c("variables", "post.match.smd", "abs_effsize")) %>%
  pivot_longer(-variables, names_to = "name", values_to = "value") %>%
  mutate(name = ifelse(name == "post.match.smd", "Matched", "All"))

order <- data2 %>%
  filter(name == "Matched") %>%
  arrange(desc(value)) %>%
  pull(variables)

data2 <- data2 %>%
  mutate(sign_positive = ifelse(
    variables == "Service Organization/State Agency" & name == "All",
    1,
    0
  ))

ggplot(data2) + geom_point(aes(
  x = variables %>% factor(levels = rev(order)),
  y = value,
  shape = name,
  color = factor(sign_positive)
)) +
  scale_shape_manual("", values = c(1, 18)) + scale_colour_manual("", values = c("black", "red"), guide =
                                                                    "none") +
  #geom_hline(yintercept = 0.05, alpha = 0.5) +
  xlab("Covariate") + ylab("Absolute Standardized Mean Difference") +
  theme(
    plot.title = element_text(hjust = 0.5),
    panel.border = element_rect(
      colour = "black",
      fill = NA,
      size = 0.75
    ),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = "black")
  ) + labs(fill = 'NEW LEGEND TITLE', color = '') +
  coord_flip()

ggsave(
  "./Replication Files/Figures/Figure 8.png",
  device = "png",
  width = 7,
  height = 3.5
)

################################################################################
### FIGURE 9: FRACTION OF APPEALS ADVANCED ON DOCKET WITHOUT DOCUMENTED CAUSE ##
################################################################################

matched_summary_aod <- matched_dataset %>%
  filter(no_aod_criteria == 1 & nodyear < 2013) %>%
  group_by(nodyear, CI) %>%
  summarise(
    mean_aod = mean(aod),
    sd = sd(aod),
    Observations = n()
  )

ggplot(data = matched_summary_aod, aes(
  x = nodyear,
  y = mean_aod,
  group = CI,
  colour = factor(CI)
)) +
  geom_point(aes(size = Observations)) + geom_line() +
  theme(
    plot.title = element_text(hjust = 0.5),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = "black")
  ) +
  labs(
    title = "Fraction of Appeals Advanced without Documented Cause",
    x = "Notice of Disagreement Year",
    y = "Fraction of Appeals Advanced",
    color = " "
  ) +
  scale_color_discrete(labels = c("No Intervention", "Congressional Intervention")) +
  scale_x_continuous(breaks = seq(2004, 2012, 4))

ggsave(
  "./Replication Files/Figures/Figure 9.png",
  device = "png",
  width = 8,
  height = 5
)

################################################################################
###### FIGURE 10: KAPLAN-MEIER CURVE FOR ADJUSTED TIME TO FINAL DECISION #######
################################################################################

fit <- survfit(Surv(adj_timeline, decided) ~ CI, data = matched_dataset)

ggsurvplot(
  fit,
  data = matched_dataset,
  conf.int = TRUE,
  legend.labs = c("No Intervention", "Congressional Intervention"),
  legend.title = "",
  xlab = "Time (Days)",
  xlim = c(0, 1440),
  break.x.by = 180,
  ylab = "Probability No Decision",
  title = "Kaplan-Meier Curve - Post-Intervention Time to Decision\n",
  font.title = c(12)
)

ggsave(
  "./Replication Files/Figures/Figure 10.png",
  device = "png",
  width = 7.5,
  height = 5
)

################################################################################
######################### FIGURE 12: MEDIATION ANALYSIS ########################
################################################################################

matched_dataset$aod_criteria <- 1 - matched_dataset$no_aod_criteria

# MEDIATION: AOD CRITERIA DOCUMENTATION
mediator_model_a <- lm(aod_criteria ~ CI, matched_dataset)
outcome_model_a <- lm(aod ~ CI + aod_criteria, matched_dataset)
results_a <-
  mediate(
    mediator_model_a,
    outcome_model_a,
    treat = 'CI',
    mediator = 'aod_criteria' ,
    boot = TRUE,
    sims = 500
  )
summary(results_a)
plot(results_a)

sens.out_a <-
  medsens(results_a,
          rho.by = 0.1,
          effect.type = "indirect",
          sims = 100)

# MEDIATION : ISSUE DOCUMENTATION
mediator_model_b <- lm(n_issues ~ CI, data = matched_dataset)
outcome_model_b <-
  lm(n_issues_allowed ~ CI + n_issues, data = matched_dataset)
results_b <-
  mediate(
    mediator_model_b,
    outcome_model_b,
    treat = 'CI',
    mediator = 'n_issues' ,
    boot = TRUE,
    sims = 500
  )
summary(results_b)
plot(results_b)

sens.out_b <-
  medsens(results_b,
          rho.by = 0.1,
          effect.type = "indirect",
          sims = 100)

png(
  file = "./Replication Files/Figures/Figure 12.png",
  width = 9,
  height = 2.5)

par(mfrow = c(1, 2), mar = c(0,0,0,0), mai = c(1, 1.5, 0.4, 0.1), tck = -0.05)
plot(
  results_a,
  main = "AOD Criteria",
  labels = c("Indirect Effect ", "Direct Effect ", "Total Effect ")
)
plot(
  results_b,
  main = "Number of Issues",
  xlim = c(0, 0.1),
  labels = c("Indirect Effect ", "Direct Effect ", "Total Effect ")
)
abline(v = c(0),
       lty = 2,
       col = "black")
par(mfrow = c(1, 1))
dev.off()

################################################################################
#####################  EFFECT OF COMMITTEE MEMBERSHIP ##########################
################################################################################

matched_dataset <-
  matched_dataset %>%
  mutate(
    congress_first_corr = case_when(
      first_corr_date >= as.Date("2003-01-07") &
        first_corr_date <= as.Date("2005-01-03") ~ 108,
      first_corr_date >= as.Date("2005-01-04") &
        first_corr_date <= as.Date("2007-01-03") ~ 109,
      first_corr_date >= as.Date("2007-01-04") &
        first_corr_date <= as.Date("2009-01-05") ~ 110,
      first_corr_date >= as.Date("2009-01-06") &
        first_corr_date <= as.Date("2011-01-04") ~ 111,
      first_corr_date >= as.Date("2011-01-05") &
        first_corr_date <= as.Date("2013-01-02") ~ 112,
      first_corr_date >= as.Date("2013-01-03") &
        first_corr_date <= as.Date("2015-01-05") ~ 113,
      first_corr_date >= as.Date("2015-01-06") &
        first_corr_date <= as.Date("2017-01-02") ~ 114,
      first_corr_date > as.Date("2017-01-03") ~ 115,
      .default = NA_integer_
    )
  )

matched_dataset <- matched_dataset %>%
  mutate(cd_match = str_pad(
    as.character(cd_match),
    width = 4,
    # match format in Census ZIP-CD crosswalk
    side = "left",
    pad = "0"
  ))

committee_info <-
  select(district_info,
         c("congress",
           "cd_match",
           "house_representative",
           "VA_committee"
         ))
committee_info <- committee_info %>%
  mutate(congress_first_corr = congress)

matched_dataset_district_info <-
  left_join(matched_dataset,
            committee_info,
            by = c("congress_first_corr", "cd_match"))
matched_dataset_district_info <-
  matched_dataset_district_info %>% mutate(VA_committee = ifelse(VA_committee ==
                                                                   2, 1, VA_committee))

t_test(
  matched_dataset_district_info %>% filter(CI == 1),
  adj_timeline ~ VA_committee,
  detailed = TRUE
)
t_test(matched_dataset_district_info,
       adj_timeline ~ VA_committee,
       detailed = TRUE)

t_test(matched_dataset_district_info %>% filter(CI == 1),
       aod ~ VA_committee,
       detailed = TRUE)
t_test(matched_dataset_district_info %>% filter(CI == 0),
       aod ~ VA_committee,
       detailed = TRUE)
t_test(matched_dataset_district_info, aod ~ VA_committee, detailed = TRUE)
